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Abstract 



Numerical simulations of the quarkonium spin splittings are done in the 
framework of lattice nonrelativistic quantum chromodynamics (NRQCD). At 
leading order in the velocity expansion the spin splittings are of 0(Mqv ), 
where Mq is the renormalized quark mass and v 2 is the mean squared quark 
velocity (v^ ~ .3 and v\ ~ .1). A systematic analysis is done of all next- 
to-leading order corrections. This includes the addition of 0{Mqv & ) rela- 
tivistic interactions, and the removal of 0{a 2 MQV i ) discretization errors in 
the leading-order interactions. Simulations are done for both S- and P-wave 
mesons, with a variety of heavy quark actions and over a wide range of lattice 
spacings. Two prescriptions for the tadpole improvement of the action are 
also studied in detail: one using the measured value of the average plaquette, 
the other using the mean link measured in Landau gauge. Next-todeading or- 
der interactions result in a very large reduction in the charmonium splittings, 
down by about 60% from their values at leading order. There are further 
indications that the velocity expansion may be poorly convergent for charmo- 
nium. Prelimary results show a small correction to the hyperfine splitting in 
the Upsilon system. 
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I. INTRODUCTION 



Quarkonium physics has been the subject of renewed theoretical interest in recent years. 
The rich phenomenology of the charmonium and Upsilon families has spurred the develop- 
ment of nonrelativistic quantum chromodynamics (NRQCD), an effective field theory that 
relies on an expansion of the action in the mean squared velocity v 2 of the heavy quarks 
(vl ~ .3 and v\ ps .1). 

NRQCD has been formulated both on the lattice and in the continuum M . Lattice 
simulations of the Upsilon and charmonium systems have recently been done by the NRQCD 
collaboration f§-|3[]. Results have also been reported for heavy-light || and be spectra [|7|], 
and some unquenched simulations have also been done ||. 

There are two key theoretical ingredients underlying lattice NRQCD calculations. One 
is an expansion of the effective action to a sufficient order in the heavy quark velocity, so 
as to obtain results of the desired accuracy. The other key ingredient is the use of tadpole 
renormalization JTIJ of the operators in the lattice action, which may then allow for reliable 
tree-level matching of the lattice theory with continuum QCD. A related problem is the 
development of sufficiently accurate discretizations of the relevant operators. 

The quarkonium spin structure is particularly sensitive to the details of the NRQCD 
Hamiltonian. The Upsilon and charmonium || spin splittings were recently analyzed 
to leaving order in v 2 by the NRQCD collaboration. At leading order in the velocity ex- 
pansion the spin splittings are of 0(Mgt> 4 ), where Mq is the renormalized quark mass. 
In the charmonium system however appreciable next-to-leading-order effects are expected, 
given the large mean squared velocity. Indeed a recent analysis of the charmonium hyper- 
fine splitting using the relativistic Fermilab action gives a result (sa 70 MeV) |l !| that is 
significantly smaller than is obtained from leading-order NRQCD (« 96 MeV) ||. 

In this paper a systematic analysis is done of all next-to-leading order corrections to 
the spin splittings. This includes the addition of 0(Mqv 6 ) relativistic interactions, as well 
as the removal of 0(a 2 MQV A ) discretization errors that are present in the leading-order 
spin-dependent operators considered in Refs. 

Results from simulations with a variety of heavy quark actions and over a wide range 
of lattice spacings are presented for the charmonium S'-wave hyperfine splitting. Some 
preliminary results are also reported for the charmonium P-wave fine structure, and for the 
Upsilon hyperfine splitting. 

Furthermore two prescriptions for defining the tadpole improvement of the action are 
studied in detail: one using the measured value of the average plaquette, as considered in 
Refs. HI!, the other using the mean link measured in Landau gauge. Landau gauge tadpole 
improvement has recently been shown to yield smaller discretization errors in the gluonic 
action (as measured by violations of rotational invariance in the heavy quark potential), 
compared to calculations using the average plaquette as input JT2]]. It is also interesting to 



note that the mean link is maximized in Landau gauge, so this prescription provides a lower 
bound on the tadpole renormalization compared to mean link determinations using other 
gauge fixings. 

An important aspect of these simulations is the removal of leading discretization errors 
in the gluonic action as well as in the NRQCD action. Specifically, an 0(a 4 )-accurate 
gluonic action is used, together with 0(a 4 )-accurate clover fields and covariant derivatives 
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in the heavy quark action. Tadpole improvement of both the gluonic and the heavy quark 
actions has recently been shown to give a good description of the spin-averaged charmonium 



spectrum even on coarse lattices with spacings a as large as .4 fm In this paper the 
charmonium spin splittings are computed on lattices with spacings in the range of about 
.17 fm to .39 fm. 

The next-to-leading order interactions are shown to result in a very large reduction in the 
charmonium hyperfine splitting, down by about 60% from the leading order result reported 
in Ref. on a lattice of comparable spacing (when the same tadpole improvement scheme 
is used); results for the triplet P-wave meson masses show that the fine structure splittings 
are also reduced by about 60%. The next-to-leading order charmonium hyperfme splitting is 
about (55±5) MeV, compared to the relativistic Fermilab action result of about (70±3) MeV 



TTJ. While the two calculations have different systematic errors, this comparsion suggests 
that further relativistic corrections, beyond the next-to-leading order considered here, are 
again very large. 

These results indicate that the NRQCD velocity expansion may be poorly convergent for 
charmonium, with the first three terms in the expansion for the hyperfme splitting apparently 
oscillating in sign. Another possibility is that there are large radiative corrections to the 
coefficients of the spin-dependent operators in the effective action. Preliminary results for 
the Upsilon hyperfine splitting at next-to-leading order in v 2 (and with leading discretization 
errors removed) show little change from the leading-order result reported in Ref. [7j . 

The dependence of the hyperfine splitting on lattice spacing is also analysed. The split- 
ting shows very little dependence on a when Landau gauge tadpole renormalization is used. 
The average plaquette tadpole scheme on the other hand has large discretization errors, but 
the results are not inconsistent with an extrapolation to the same splitting at zero lattice 
spacing in the two schemes. These results provide further support for the use of Landau 



gauge tadpole renormalization [12 



II. QUARK AND GAUGE-FIELD ACTIONS 

The NRQCD effective action is based on power-counting rules for the magnitude of heavy 
quark and gauge field operators in quarkonium states. The expansion parameters are the 
mean squared velocity of the heavy quarks in the bound state, and the strong coupling 
constant. The coefficients of the operators in the effective action, to a given order in v 2 , are 
determined by matching the predictions of NRQCD with those of full QCD |1]||. 

The heavy quark lattice Hamiltonian is conveniently decomposed into the leading covari- 
ant kinetic energy operator H , plus relativistic and discretization corrections 5H. Following 
Refs. |5|,[|, the quark Green function is given by 

G t+1 =(l-^)V|(l-^) n (l-a5i/)G, (t>l), (1) 
where the initial evolution is set by 

The kinetic energy operator H is of 0(v 2 ), and is given by 
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where M c ° is the bare quark mass and A^ is the lattice Laplacian. Relativistic corrections 
are organized in powers of the heavy quark velocity, with terms up to 0(v 6 ) considered here: 



5H = 5H {4) + 5H {6) . (4) 
Simulations with the complete set of 0(t> 4 ) corrections were reported in Refs. @.[|: 




terms come from finite lattice spacing corrections to the lattice Laplacian and the lattice 
time derivative respectively The parameter n is introduced to remove instabilities in the 
heavy quark propagator caused by the highest momentum modes of the theory ||. 

While 6H^ yields spin-averaged spectra to next-to-leading order in the velocity expan- 
sion, it contains only leading order spin-dependent interactions. In potential model language 
the third term above (03) generates the spin-orbit and tensor potentials which drive the P- 
wave fine structure, while the fourth term (04) generates the color-magnetic S'-wave hyperfine 
splitting. 

Spin-dependent interactions at 0(t> 6 ) were derived in Ref. 0: 



Notice the field strength bilinear (eg), which is peculiar to the nonAbelian theory. There are 
additional 0(v 6 ) terms which contribute to spin-averaged spectra; these are not considered 
here. 

The derivative operators and the fields are evaluated with their leading discretization 
errors removed, in order to minimize the effects of lattice artifacts on the spin splittings. 
This is indicated by the tilda superscripts on these operators in Eqs. (||) and (||). At leading 
order in a the action of the symmetric lattice derivative Aj is defined by 




-c 9 



8(M°) 



-(7-ExE. 



(6) 



aAiG(x) = \[Ui{x)G{x + ai) - U}(. 



x — ai)G(x — ai)]. 



(7) 



At tree level the 0(a 4 )-accurate derivative operator Aj that is used in Eqs. fl5|) and (||) is 
given by @ 



a^a.-^a.Waa^, (8) 

o 



where A^ and A - "* are leading-order forward and backward covariant finite differences: 

aA^G^x) = Ui(x)G(x + ai) - G(x), 

a/sf ) G(x) = G(x) - Uj(x - az)G{x - ai). (9) 

At leading-order in a the lattice Laplacian A^ 2 ) is expressed in terms of covariant second 
order differences 

A( 2 )=£Af\ (10 ) 

i 

where 

a 2 Af ] G(x) = Ui(x)G(x + ai) - 2G(x) + Uj(x - ai)G(x - at). (11) 
The 0(a 4 )-accurate Laplacian A^ 2 ) is used in Eq. (|]), and at tree- level is given by 

A( 2 ) = A( 2 ) - -A< 4 \ (12) 

12 K J 

where 

A (4) =E(Ai 2) ) 2 (13) 



is a lattice representation of the continuum operator Yli Df. Note that discretization errors 
in Hq have been removed in this way by the addition of the C5 term in Eq. (|5|). What is 
new in the present work is the removal of discretization errors in the leading order P-wave 
interaction (03), and the Darwin term (02), as well as the addition of the 0(v 6 ) spin terms 
in Eq. (H). 

To complete the removal of discretization errors from the spin-dependent interactions 
0(a 4 )-accurate chromo-electric and -magnetic fields have been used in these simulations. In 
Refs. the leading order field strength F^ was evaluated using the standard traceless 
clover operator: 

F, u {x) = i (Sl^x) - nl(x)) - ~Im(Tr£V(z)) , (14) 

where fl^ u is an average over the 4 counterclockwise plaquettes in the (/x, v) plane containing 
the site x 

n ^ = -7 H E U a (x)Up(x + a)U- a (x + a + p)U^{x + p). (15) 

4 a=±n 0=±v 

The 0(a 4 )-accurate field strength F^ v used here is computed following the analysis of Ref. 
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5 1 

U fl (x)F lxu (x + fyWJx) 



F^ v (x) ^F^ v (x) 



+Ufa - fi)F^(x - fi)U^{x -p)-(p.++ i/)] . (16) 

The "improved" chromo-electric and -magnetic fields are defined in terms of the field 
strength, E { = F 4i and B { = \e ijk F jk . 

It should be noted that some of the discretization errors that are removed by using A, 
hS 2 \ and Ffj, u in Eqs. (|) and (|), in place of their leading order counterparts, may in fact 
be comparable to or smaller than higher order relativistic corrections that are not included 
here. This includes for example the use of A® and F^ in the O (m 6 ) terms in Eq. @, and 
the use of the improved operators in the Darwin term (02). On the other hand, the use of 
improved operators in the leading spin-dependent interactions (03 and C4) corrects for errors 
of 0(a 2 M c v 4 ) in the spin splittings which, for the range of lattice spacings studied here, may 
be comparable to the 0(M c v 6 ) contributions from 5H^ 6 \ 

At tree-level all of the coefficients Cj in Eqs. (|5|) and (|6]) are one. However, very large 
radiative corrections in the lattice theory can arise from tadpoles that are induced by the 
nonlinear connection between the link variables and the continuum gauge fields. Most 



of the effects of tadpoles can be removed by a mean-field renormalization of the link flTOf: 



U,{x) - (17) 

Uq 

The links are rescaled in the simulation before they are input to the quark propagator 
subroutine, to be sure that Eq. (|T7|) is correctly implemented in all terms in the heavy 
quark action. 

In most previous work the fourth root of the average plaquette has been used to set the 
value of Mo: 

u 0tP = (iReTrf/ pl ) 1/4 . (18) 

Simulations were done here with this renormalization prescription. In addition, simulations 
were also done using the mean link in Landau gauge to set u , as recently suggested by 



Lepage [12 



M 0)i = (|ReTrC/ M ), <W = 0, (19) 



where a standard lattice implementation of the continuum Landau gauge fixing is used 
|T4| (it was found that the removal of leading discretization errors in the lattice version of 
<9 M A M = results in a negligible change to the value of mo,l)- 

Finally, the gauge-field configurations were generated using an 0(a 4 )-accurate tadpole- 



improved action [13| 



S[U] = |ReTr (l - U pl ) - £ |ReTr (1 - U lt ) , (20) 
pi rt 

where the sums are over all oriented lxl plaquettes and 1x2 rectangles. 
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III. MESON PROPAGATORS 



In order to increase the overlap of the meson propagators with the ground states of 



interest here, a gauge- covariant smearing procedure has been used [15|. A meson creation 



operator is constructed from quark and antiquark creation operators ip^ and 0)11)11 : 

5>t (f)r(f)x t (f); ( 2i) 

X 

with 

r(£)=fi(£) 7 (£), (22) 

where Q(x) is a 2 x 2 matrix in spin-space, with derivative operator entries, which gives the 
quantum numbers of the state of interest. 7(2?) is a gauge-covariant local smearing operator, 
which is taken to have the simple form |15] 



1 (x) = (l + eA^(x)) ns (23) 

(an invariant under the lattice cubic group). The weight e and the number of smearing 
iterations n s are adjusted to optimize the overlap with the ground state. 

The meson correlation function G meS on at zero momentum is then given by 

G mcson (p = 0,t) = ]TTr [G\{y - x)T\ sk) (g}G t (y - x)T {sc) (x)} , (24) 
y 

where different smearing parameters may be used at the source and sink, and where a single 
spatial origin x for the meson propagator was generally used. Finite momentum propagators 
for the 1 S'o were analyzed using a local 5-function source and sink: 

G mcson (p,t) = X> [G\{y- x)G t (y- x)} e ~^~ s \ (25) 
y 

Correlation functions were computed for the 1 S Q (£l = I), 3 Si (tt = a t ) and 1 P 1 (VL = A { ) 
mesons. The three triplet P-wave correlators ( 3 Po) 3 -Pb 3 -p2) were also analyzed; the relevant 
operators Q for these states are tabulated in Ref. ||. Only selected meson polarizations were 
used. Propagators were generated for all (equal) quark-antiquark colors but, in order to save 
computer time, the initial quark and antiquark spins were set to 1. The 3 Siz and 1 S'o states 
were thus obtained from a single propagator, since the lxl component of the spin matrix 
Q is the same for both states 0. All three polarizations of the 1 Pi were generated, but only 
one from each of the 3 Pi and 3 P2- 



IV. RESULTS 



Three lattices were generated using the mean link in Landau gauge to set the tadpole 
factor (mo,l) and four lattices with comparable spacings were generated using the average 
plaquette tadpole («o,p) • The parameters of these seven lattices are given in Tables p] and 0. 
In order to distinguish between the two sets of simulations, (3l is used to denote the lattice 
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coupling when Landau-gauge tadpole renormalization is used, and (3p when the average 
plaquette is used. 

A standard Cabbibo-Marinari pseudo-heat bath was used to generated the gauge field 
configurations. Integrated autocorrelation times r- irit were checked for all correlation func- 
tions, and were found to be remarkably short; 10 updates between measurements yields 
Tint ~ 0.5 on the three lattices with a < .2 fm, and on the coarser lattices 5 updates was 
found to be sufficient. 

Smeared-smeared correlators were used for the P-waves, while local sources and smeared 
sinks were used for the S-waves. Ten smearing iterations [n s = 10 in Eq. fl2"3"|)] were used 
for the three lattices with a < .2 fm, 5 iterations for the lattices with a near .28 fm, and 
2 iterations for the lattices with a ~ .39 fm. A smearing weight e = 1/12 was used in all 
cases. 

The lattice spacings are determined from the spin- averaged IP — IS mass difference, 
following Refs. This mass difference is known to be independent of the quark mass in 

the charm to bottom region. For this purpose the singlet 1 P\ and the spin-averaged 3 Si, 1 S 
masses were used. The simulation results for the splitting were fixed to the experimental 
value for charmonium of 458 MeV. 

After the lattice spacing was extracted, the kinetic mass M kin (in physical units) of the 
1 S state was determined by fitting the energy P P of the boosted state to the form 

p2 

E P - E = . (26) 

Fits were made to the state with momentum components (1,0,0) in units of 2ir/(Na); 
in some cases simultaneous fits including states with momentum components (1, 1, 0) and 
(1, 1, 1) were also done, with little change to the fit values of Mk m . A dispersion relation 
including relativistic corrections was also tried, and the resulting changes to the fit 
values of M kin were within a few percent, as expected on these lattices. 

The correct values of the bare charm mass M c ° were determined by tuning so that M kin 
agrees with the experimental value of the mass of the i] c (2.98 GeV). The bare masses are 
listed in Tables | and [T| and all yield M kin = 2.9(1) GeV. For bare masses aM? < 1.2 
a stability parameter n = 4 was used in the quark propagators, Eqs. (|l|) and (|||); for 
1.2 < aM c ° < 1.5 n = 3 was used, and for the larger bare masses n = 2 was used. 

Effective mass plots m e g(T) = — log(G(T)/G(T — 1)) for several lattices are shown in 
Figures using jackknife errors. Single exponential fits to the correlation functions are 
used to get the best estimates of the masses of the individual states. The fitting procedure 
included the full covariance matrix for the data, using the svd algorithm [|l|| . The correlation 
functions for states of a given partial wave are highly correlated; following Refs. |5||| a spin 
splitting 5E was obtained from a correlated fit of the form 

^ meson, A (^) , 

G mcsoa , B (t)=c B e~ {EA+5E)t . (27) 

Detailed fit results for several lattices are reported in Tables [1TT|-[VT|. The statistical 
errors were estimated using bootstrap ensembles of 1000 samples. Final estimates of the 
dimensionless energies are obtained from these fits by finding two or three successive t m in/^max 
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intervals for which the fit results overlap within statistical errors; acceptable Q values were 
obtained in all cases at these t m i n /t max values. Estimates of the systematic errors in the final 
fit results are taken from the largest statistical errors in the overlapping intervals. 

The final fit results are shown in Tables |VII| and [VIII| , where the resulting lattice spacings 
and hyperfine splittings in physical units are also given. The dominant error in the splitting 
comes from the systematic error in the determination of the bare quark mass. The error in 
the mass comes in part from the uncertainty in a -1 , which has been included in the error 
estimates for the splittings in physical units; however, there is a further systematic error of 
order 10% in the quark mass determination, coming from higher order (spin-independent) 
relativistic corrections |E||. 



V. DISCUSSION 

The hyperfine splittings are plotted as a function of lattice spacing in Fig. |5|, where the 



results from the relativistic Fermilab action [11] and the leading order NRQCD calculation 



U are included. Some coarse lattice results from the tadpole-improved relativistic D234 [T7 
action are also shown. 

The next-to-leading order corrections result in a very large reduction in the hyperfine 
splitting, down by about 60% from the leading order result on a lattice of comparable 
spacing, when the same plaquette tadpole renormalization scheme (mo.p) is used in both 
cases. 

The hyperfine splitting shows very little a dependence when the Landau gauge tadpole 
scheme «o,l is used. The results with it 0i p on the other hand have large discretization errors, 
which prevents a reliable extrapolation to zero lattice spacing in this case; however, the 
results are not inconsistent with an extrapolation to the same splitting as is obtained with 
Landau gauge mo,l- 

From these results the hyperfine splitting at next-to-leading order in the velocity ex- 
pansion, and at zero lattice spacing, can be estimated at roughly (55 ± 5) MeV. This can 



be compared with the Fermilab action result of approximately (70 ± 3) MeV JTTJj. While 
the two calculations have different systematic errors, this suggests that further relativistic 
corrections, beyond the next-to-leading-order considered here, are again very large. [The 
experimental value is (118 ± 2) MeV, which indicates that there are significant effects due 



to quenching \\ 1 



These considerations are supported by results obtained here for the triplet P-wave spectra 
( 3 Po, 3 -Pi) 3 -£2)- The next-to- leading order fine structure splittings are apparently very small, 
and much better statistics are required for an accurate determination. Simulations with mo,l 
give a 3 P 2 — 3 Po splitting of about (30 ±15) MeV, down by about 60% from the leading-order 
result of (110 ± 10) MeV reported in Ref. || (the experimental value is (141 ± 10) MeV). 
The next-to-leading order fine structure splittings with u Q> p are even harder to measure, but 
the results suggest that the splittings may actually be in the wrong order with that tadpole 
scheme (as least in the range of lattice spacings analyzed here). 

These results indicate that the NRQCD velocity expansion for charmonium may be 
poorly convergent, with the first three terms in the expansion for the hyperfine splitting 
apparently oscillating in sign. However there are other sources of systematic error in the 
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NRQCD action which must also be considered. Radiative corrections to the operator coef- 
ficients C3 and C4 in Eq. @ are of particular importance. 

It is worthwhile to assess the relative importance of the various next-to-leading order 
corrections that have been considered here. Within a given tadpole renormalization scheme 
the most important correction for the charmonium system comes from the relativistic spin- 
dependent interactions SH^ [Eq. (|6])], and these drive the large reduction in the splittings. 
There is some indication that the field strength bilinear (term eg) plays a relatively small 
role in these effects. The use of 0(a 4 )-accurate clover fields increases the spin splittings, a 
correction amounting to about 20% of the 0(a 2 )-accurate splittings on the coarsest lattices 
considered here (this correction falls below about 10% at the smallest spacings). 

The effect of a change in the tadpole renormalization scheme is very significant. The spin 
splittings vary as 1/mq due to the renormalization of the clover field interactions (except for 
the field strength bilinear, which varies as 1/mq). This renormalization causes most of the 
change due to tadpole scheme, as can be seen from the values of m ,l and w 0i p i n Tables | 
and H (additional changes in the splittings are presumably caused by the renormalization 
of the gluonic action). For example, at a ~ .18 fm the ratio (wo,p/^o,l) 4 is about 1.2, and 
at a ~ .39 fm the ratio is about 1.5. 

Relativistic corrections are expected to be much smaller for the Upsilon system. Prelim- 
inary results from a next-to- leading order calculation at (3p = 7.2 give an T — rjf, hyperfine 
splitting of (22.4 ± 1.3) MeV, using a bare mass aM' = 3.15. This is within errors of the 
leading-order result reported in Ref. on a lattice with comparable spacing (using wo,p), 
which is consistent with the velocity expansion, given the fact that v\ ~ 0.1. The small size 
of the net correction, relative to the charmonium system, also appears to be driven by a near 
cancellation of the 0(a 2 M<gt> 4 ) discretization corrections (which tend to raise the splittings) 
and the 0(Mqv 6 ) relativistic corrections (which tend to lower them), which are much closer 
in magnitude in the Upsilon system. This was demontrated by doing a calculation with a 
Wilson gauge field action at (3 = 5.7, without removing discretization errors in the clover 
field, but including the 0(Mqv & ) relativistic interactions; this results in a reduction of the 
Upsilon hyperfine splitting by about 15% compared to the leading order calculation @ at 
the same (3. 

It is also interesting to note that the lattice spacing as determined from the IP — 15 
splitting is different for Upsilon and charmonium. At (3p = 7.2 preliminary results give 
ax = 0.146(9) fm, compared to = 0.171(4) fm. This is comparable to the difference 
between the two determinations of the spacing in Ref. |7J . 

VI. SUMMARY AND OUTLOOK 

It has been shown that spin-dependent interactions at next-to-leading order in the 
NRQCD velocity expansion yield very large corrections to the charmonium spin splittings, 
down by about 60% from their values at leading order (when the same tadpole-improvement 
scheme is used on lattices with comparable spacings). There are indications that further 
relativistic corrections for charmonium are also very large. The corrections to the hyperfine 
splitting in the Upsilon system are small. More work needs to be done in order to assess 
the validity of the NRQCD effective action in simulations of charmonium, including better 
measurements of the triplet P-wave spectra. Estimates of the radiative corrections to the 
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operator coefficients in the heavy quark action are also needed. More complete calculations 
of the Upsilon splittings would also provide useful information. The results obtained here 
provide further support for the use of Landau gauge tadpole renormalization. 
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TABLES 







(iReTrC/ pl )V4 


a (fm) 


oM c ° 


Volume 


-Vmcas 


7.4 


.829 


.875 


.18 


1.18 


10 3 x 16 


752 


7.0 


.780 


.850 


.28 


1.90 


6 3 x 10 


1740 


6.6 


.743 


.825 


.39 


2.65 


6 3 x 10 


2092 



TABLE I. Simulation parameters using the Landau gauge mean link to determine the tadpole 
renormalization, uo,l = (jReTrf/^) (second column). -/V mcas is the number of configurations used 
for propagator measurements. 



0P 


(iReTV*7 M ) 


(iReTr[/ pl )V4 


a (fm) 


aM c ° 


Volume 


-Vmeas 


7.2 


.834 


.874 


.17 


0.81 


10 3 x 16 


474 


7.0 


.810 


.865 


.21 


1.10 


8 3 x 10 


923 


6.8 


.786 


.854 


.26 


1.43 


6 3 x 10 


1815 


6.25 


.738 


.821 


.39 


2.30 


6 3 x 10 


2841 



TABLE II. Simulation parameters using the average plaquette to determine the tadpole renor- 
malization, uq,p = (|ReTr t/'p^) 1 / 4 (third column). 



^min/^max 


'Pi 


3 s, 


1 S'o 


3 c 1 Q 
Di — Do 


2/16 


0.735(5) 


0.298(2) 


0.242(1) 


0.0537(6) 


3/16 


0.715(6) 


0.293(2) 


0.239(1) 


0.0533(5) 


4/16 


0.709(9) 


0.289(2) 


0.237(1) 


0.0518(6) 


5/16 


0.694(12) 


0.286(2) 


0.235(1) 


0.0506(6) 


6/16 


0.688(16) 


0.285(2) 


0.235(1) 


0.0497(7) 


7/16 


0.672(23) 


0.284(2) 


0.234(1) 


0.0500(8) 


8/16 


0.661(30) 


0.284(2) 


0.234(1) 


0.0498(9) 


9/16 


0.700(51) 


0.284(2) 


0.234(1) 


0.0494(10) 



TABLE III. Examples of fits to Landau-gauge tadpole simulation at /3l = 7.4 (a = .18 fm). 
Single exponential fits were used for each individual state, and a correlated 8E fit for the 3 5i and 
1 S'o was used to get the hyperfine splitting. 
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^min/^m: 

2/10 
3/10 
4/10 
5/10 
6/10 
7/10 



1.280(8) 

1.268(15) 

1.279(35) 

1.198(71) 

1.18(14) 

1.44(43) 



3 ^i 
0.404(1) 
0.395(1) 
0.393(1) 
0.392(2) 
0.392(2) 
0.392(2) 



1 '5'o 
0.304(1) 
0.299(1) 
0.298(1) 
0.298(1) 
0.298(1) 
0.298(1) 



3 c 1 c 

0.0980(6) 

0.0969(6) 

0.0950(7) 

0.0945(7) 

0.0946(10) 

0.0944(12) 



TABLE IV. Examples of fits to Landau-gauge tadpole simulation at /3l = 6.6 



Fits were done as in Table III. 



.39 fm). 



^minAmax 1 J~\ jSl jSq 3 S\ — 1 Sp 

2/16 1.120(6) 0.715(2) 0.675(2) 0.0399(5) 

3/16 1.101(7) 0.708(2) 0.670(2) 0.0395(5) 

4/16 1.089(9) 0.706(2) 0.669(2) 0.0369(5) 

5/16 1.083(12) 0.704(2) 0.667(2) 0.0363(6) 

6/16 1.087(15) 0.702(2) 0.667(2) 0.0357(7) 

7/16 1.083(21) 0.700(2) 0.666(2) 0.0352(7) 

8/16 1.093(29) 0.701(2) 0.666(2) 0.0352(8) 

TABLE V. Examples of fits to average plaquette tadpole simulation at (3p = 7.2 (a = .17 fm). 



Fits were done as in Table [IL 



'mill / ' ma 



2/10 
3/10 
4/10 
5/10 
6/10 
7/10 



1.773(7) 
1.754(14) 
1.754(26) 
1.771(60) 



3 5i 
0.858(1) 
0.853(1) 
0.852(1) 
0.851(1) 
0.851(1) 
0.850(1) 



0.801(1) 
0.797(1) 
0.795(1) 
0.795(1) 
0.795(1) 
0.795(1) 



3 c 1 c 
£1 — <->0 

0.0575(3) 

0.0570(3) 

0.0565(3) 

0.0561(4) 

0.0560(4) 

0.0556(6) 



TABLE VI. Examples of fits to average plaquette tadpole simulation at (3p = 6.25 (a = .39 fm). 



Fits were done as in Table til. 
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Pl 


'Pi 


'S 


3 5i - % 


a 


Hyperfine 










(fm) 


(MeV) 


7.4 


0.68(2) 


0.234(1) 


0.0497(7) 


0.176(9) 


55.7(27) 


7.0 


1.01(2) 


0.305(1) 


0.0772(8) 


0.278(9) 


54.6(18) 


6.6 


1.27(3) 


0.298(1) 


0.0945(7) 


0.388(13) 


48.0(17) 



TABLE VII. Final fit results for the dimensionless energies from Landau-gauge tadpole simu- 
lations; the resulting lattice spacings and hyperfine splittings in physical units are also shown. The 
quoted errors in the hyperfine splittings in physical units include the systematic errors in a. 



Pp 


'Pi 


'So 


3 S 1 - 'So 


a 


Hyperfine 










(fm) 


(MeV) 


7.2 


1.09(1) 


0.666(2) 


0.0352(7) 


0.171(4) 


40.5(13) 


7.0 


1.23(2) 


0.728(2) 


0.0354(7) 


0.205(9) 


34.1(16) 


6.8 


1.42(2) 


0.790(2) 


0.0427(6) 


0.257(9) 


32.7(12) 


6.25 


1.75(3) 


0.795(1) 


0.0560(6) 


0.393(13) 


28.1(10) 



TABLE VIII. Final fit results from average plaquette simulations. 
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FIGURES 



FIG. 1. Effective mass plot for j3 L = 7.4 (a = .18 fm): l P x state (□) and l S Q state (o). 
FIG. 2. Effective mass plot for /?£, = 6.6 (a = .39 fm): x Pi state (□) and 1 Sq state (o). 
FIG. 3. Effective mass plot for (3 P = 7.2 (a = .17 fm): 1 P 1 state (□) and state (o). 

FIG. 4. Effective mass plot for f3 P = 6.25 (a = .39 fm): 1 P 1 state (□) and l S state (o). 

FIG. 5. Hyperfine splittings versus lattice spacing squared. The next-to-leading order NRQCD 
results with Landau gauge tadpoles (■) and with average plaquette tadpoles (o) are shown, as well 
as the leading-order result (x) from Ref. M. Also shown are results obtained with the relativistic 
Fermilab action (□) [jll]], and coarse lattice results from the relativistic tadpole-improved D234 
action (A) [[l?]]. The experimental value is (118 ± 2) MeV. 



16 



ho 



ISO 



n ? 1 1 1 1 1 — *~ 



o 



o 



o 



o 



o 



o 



o 



o 



-e- 



-B- 



-0 



□ 



n 



o 



hCH 



h-CH 



I — □ — I 



I □ 1 



I □ 1 



I □ 1 



o o o o 
Kd k qd 




ho 



i\3 




II 



^ ^ — k — k — k — k — k 



to 




to 



Hyperfine (MeV) 





1 1 1 


1 1 1 
1 □ 1 


1 


• 








o 


O 


h-BH 


rAI 


• 

o 
•k 


hCH 






• 


O 






B * 




■ 












® 








• 








• 

-k 








• 

->< 


■ 

O 

1 1 1 


1 1 1 


l 



